%Written by Michael R. Springborn
%This version: Dec. 16, 2011

% Free to use and modify with proper acknowledgement. No warranties. 
%*******************************************************************************
%                         ***GENERAL DESCRIPTION***     

%This code generates maximum likelihood parameter estimates for use as starting values 
%in the Bayesian estimation implemented in runBayesEst.m
%*******************************************************************************

clear; clc
work_dir='D:\Bayesian_Analysis'
cd(work_dir)

%Select the link model for modeling binomial outcomes: 
%model='logit';
model='cauchit';

%Choose whether or not to correct for the stratified sample using Manski-Lerman weights.
MLw='no'; %Manski-Lerman weights
%MLw='yes'; %Manski-Lerman weights

%Assumed base rate of weeds (unconditional probability that Y=1).  Note: referred to in manuscript as "tau".
Q=0.05; 
%Q=0.02 

%% Load training data set
importfile('D:\data\data_WRA.csv'); %For descriptions of field names see Excel file of same name in same folder.
N=length(ID_uniqespec); %Number of observations in training data.
y=Bin_DepVar; 
X=[Score_assessor];

%% Likelihood function elements
%Manski-Lerman weights, corrective weighting for stratified sampling
ybar=mean(y);    
W1=Q/ybar;
W0=(1-Q)/(1-ybar);

%Link function: 
if strcmp(model,'logit')  %logit model
    pfun= @(X,b) 1./(1+exp(-[ones(size(X,1),1) X]*b)); %where b is 2xZ
elseif strcmp(model,'cauchit')  %Cauchit model
    pfun= @(X,b) .5 + atan([ones(size(X,1),1) X]*b)/pi; %where b is 2xZ
end
%Log-likelihood function:
if strcmp(MLw,'yes')
    lnlike = @(b,y,X,W1,W0) sum( (repmat(y,1,size(b,2))*W1 + (1-repmat(y,1,length(b(1,:))))*W0)  .*  (  repmat(y,1,length(b(1,:))).*log(pfun(X,b)) + (1-repmat(y,1,length(b(1,:)))).*log(1-pfun(X,b))   ));
elseif strcmp(MLw,'no')
    lnlike = @(b,y,X,W1,W0) sum(  repmat(y,1,size(b,2)).*log(pfun(X,b)) + (1-repmat(y,1,length(b(1,:)))).*log(1-pfun(X,b))   );
end        

%initial guess
bmle02=[-5.3 0.36]; %for Q=0.02

[B,fv,eflag] = fminsearch(@(b) -lnlike(b,y,X,W1,W0),bmle02'); 
B  %MLEs

return

%bmle:
%Logit:
%Uncorrected: bmle: .1457 .2945 
%Q=0.02, bmle=-5.2620    0.3510;
%Q=0.05, bmle=-4.2130    0.3331

%Cauchit:
%Uncorrected: bmle: 0.2190    0.3682
%Q=0.02, bmle=-31.2032  2.4559;  
%Q=0.05, bmle=-12.2617  0.9860;

%std MLE: bmle: .1457 .2945
%King and Zeng prior correction: B0-ln(stuff)  pg 1415, stat in medicine 2002.
% Q=0.02
% b0corr=.1457 - log(((1-Q)/Q)*(ybar/(1-ybar)))










